knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(terra)
library(oharac)
library(tidyverse)
library(here)
source(here('common_fxns.R'))Compare vulnerability results between unweighted means (script 2) and functional-vulnerability-weighted means (script 3).
Inputs are mean (and sd?) vulnerability maps from scripts 2 and 3.
Loop over several focal stressors. For each stressor, read in the maps of species mean \(\bar x_{spp}\) and fv-weighted mean \(\bar x_{FE}\) (and similar for coefficient of variance). Calculate the following comparisons:
squish_rast <- function(r, qtile = .999) {
r_vals <- values(r); r_vals <- r_vals[!is.na(r_vals)]
r_zlim <- max(abs(quantile(r_vals, 1 - qtile)), abs(quantile(r_vals, qtile)))
values(r)[values(r) > r_zlim] <- r_zlim
values(r)[values(r) < -r_zlim] <- -r_zlim
return(list(r = r, zlim = r_zlim))
}vuln_spp_dir <- here_anx('_output/vuln_maps/vuln_maps_by_species')
mean_fs_spp <- list.files(vuln_spp_dir, pattern = '_mean.tif', full.names = TRUE)
# sdev_fs_spp <- list.files(vuln_spp_dir, pattern = '_sdev.tif', full.names = TRUE)
vuln_fe_dir <- here_anx('_output/vuln_maps/vuln_maps_by_funct_entity')
mean_fs_fe <- list.files(vuln_fe_dir, pattern = '_mean.tif', full.names = TRUE)
# sdev_fs_fe <- list.files(vuln_fe_dir, pattern = '_sdev.tif', full.names = TRUE)
focal_strs <- c('biomass_removal', 'bycatch',
'wildlife_strike',
'habitat_loss_degradation',
'nutrient_pollution', 'ocean_acidification',
'sst_rise', 'marine_heat_wave',
'light_pollution')
for(f in focal_strs) { ### f <- focal_strs[3]
str_name <- str_replace_all(f, '[^a-z]+', ' ') %>%
str_squish() %>%
str_to_title()
message('Calculating difference maps for ', tolower(str_name), '...')
mean_f_spp <- mean_fs_spp[str_detect(basename(mean_fs_spp), f)]
# sdev_f_spp <- sdev_fs_spp[str_detect(basename(sdev_fs_spp), f)]
mean_rast_spp <- terra::rast(mean_f_spp)
# sdev_rast_spp <- raster::raster(sdev_f_spp)
# cv_rast_spp <- sdev_rast_spp / mean_rast_spp
mean_f_fe <- mean_fs_fe[str_detect(basename(mean_fs_fe), f)]
# sdev_f_fe <- sdev_fs_fe[str_detect(basename(sdev_fs_fe), f)]
mean_rast_fe <- terra::rast(mean_f_fe)
# sdev_rast_fe <- raster::raster(sdev_f_fe)
# cv_rast_fe <- sdev_rast_fe / mean_rast_fe
### Set up rasters; trim extreme values
mean_diff_rast <- (mean_rast_fe - mean_rast_spp) / mean_rast_spp
mean_diff_squished <- squish_rast(mean_diff_rast, qtile = .999)
# cv_diff_rast <- (cv_rast_fe - cv_rast_spp) / cv_rast_spp
# cv_r_sq <- squish_rast(cv_diff_rast, qtile = .999)
m_d_mask_0.10 <- m_d_mask_0.25 <- mean_diff_squished$r
values(m_d_mask_0.10)[values(mean_rast_spp) < .10 & values(mean_rast_fe) < .10] <- NA
values(m_d_mask_0.25)[values(mean_rast_spp) < .25 & values(mean_rast_fe) < .25] <- NA
message('Plotting difference maps for ', tolower(str_name), '...')
### Set up plot for diverging palette, and symmetric z limits around zero
map_cols <- hcl.colors(palette = 'Red-Green', n = 50, rev = TRUE)
# [1] "Blue-Red" "Blue-Red 2" "Blue-Red 3" "Red-Green" "Purple-Green"
# [6] "Purple-Brown" "Green-Brown" "Blue-Yellow 2" "Blue-Yellow 3" "Green-Orange"
# [11] "Cyan-Magenta" "Tropic" "Broc" "Cork" "Vik"
# [16] "Berlin" "Lisbon" "Tofino"
plot(mean_diff_squished$r, col = map_cols,
main = paste0('% Diff in mean vuln: ', str_name),
zlim = c(-mean_diff_squished$zlim, mean_diff_squished$zlim),
legend = TRUE, axes = FALSE)
plot(m_d_mask_0.10, col = map_cols,
main = paste0('% Diff in mean vuln masked 0.10: ', str_name),
zlim = c(-mean_diff_squished$zlim, mean_diff_squished$zlim),
legend = TRUE, axes = FALSE)
plot(m_d_mask_0.25, col = map_cols,
main = paste0('% Diff in mean vuln masked 0.25: ', str_name),
zlim = c(-mean_diff_squished$zlim, mean_diff_squished$zlim),
legend = TRUE, axes = FALSE)
# plot(cv_r_sq$r, col = map_cols, main = paste0('% Diff in coef of var: ', str_name),
# zlim = c(-cv_r_sq$zlim, cv_r_sq$zlim),
# legend = TRUE, axes = FALSE)
}